(in-package :maxima)

(defun maxima-repr (expr)
  (let ((string (with-output-to-string (*standard-output*)
		  (let (($display2d nil))
		    (displa expr)))))
    (subseq string 0 (1- (length string)))))

(defclass point ()
  ((x :initarg :x
      :accessor point-x)
   (y :initarg :y
      :accessor point-y)))

(defmethod print-object ((obj point) stream)
  (print-unreadable-object (obj stream :type t)
    (format stream "~A, ~A"
	    (maxima-repr (point-x obj))
	    (maxima-repr (point-y obj)))))

(defclass vec ()
  ((x-component :initarg :x-component
		:accessor x-component)
   (y-component :initarg :y-component
		:accessor y-component)))

(defun degrees->radians (degrees)
  ($ratsimp (m* degrees (m// '$%pi 180))))

(defmethod print-object ((obj vec) stream)
  (print-unreadable-object (obj stream :type t)
    (format stream "~A, ~A"
	    (maxima-repr (x-component obj))
	    (maxima-repr (y-component obj)))))

(defparameter v (make-instance 'vec
			       :x-component 3
			       :y-component 4))

(defmethod angle ((v vec))
  ($ratsimp `((%atan) ,(m// (y-component v)
			    (x-component v)))))

(defmethod magnitude ((v vec))
  ($ratsimp `((%sqrt) ,(m+ (m^ (x-component v) 2)
			   (m^ (y-component v) 2)))))

(defmethod (setf magnitude) (value (v vec))
  (let ((theta (angle v)))
    (setf (x-component v)
	  ($ratsimp (m* `((%cos simp) ,theta) value)))
    (setf (y-component v)
	  ($ratsimp (m* `((%sin simp) ,theta) value)))))

(defmethod (setf angle) (value (v vec))
  (let ((magnitude (magnitude v)))
    (setf (x-component v)
	  ($ratsimp (m* `((%cos simp) ,value) magnitude)))
    (setf (y-component v)
	  ($ratsimp (m* `((%sin simp) ,value) magnitude)))))

(defun make-vector (x y)
  "Make a vector given cartesian coordinates"
  (make-instance 'vec :x-component x :y-component y))

(defun make-vector-polar (r theta)
  "Make a vector given polar coordinates"
  (let ((x ($ratsimp (m* `((%cos simp) ,theta) r)))
	(y ($ratsimp (m* `((%sin simp) ,theta) r))))
    (make-vector x y)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defclass free-body ()
  ((acceleration :initarg :acceleration
		 :accessor acceleration)
   (velocity :initarg :velocity
	     :accessor velocity)
   (coordinates :initarg :coordinates
		:accessor coordinates
		:type point
		:initform (make-instance 'point :x 0 :y 0))
   (initial-coordinates :accessor initial-coordinates
			:initarg :initial-coordinates
			:initform (make-instance 'point :x 0 :y 0))
   (mass :initarg :mass
	 :accessor mass)
   (time :initarg :time
	 :accessor body-time
	 :initform 0)))

(defmethod set-time ((body free-body) time)
  (let ((coordinates (coordinates body))
	(xi (point-x (initial-coordinates body)))
	(yi (point-y (initial-coordinates body))))
    (setf (point-x coordinates)
	  ($ratsimp (m+ (m* time (x-component (velocity body)))
			xi)))
    (setf (point-y coordinates)
	  ($ratsimp (m+ (m* (m// 1 2)
			    (y-component (acceleration body))
			    (m^ time 2))
			(m* time (y-component (velocity body)))
			yi)))))

(defmethod set-x ((body free-body) x)
  (let ((time ($ratsimp (m// (m+ (m* -1 (point-x (initial-coordinates body)))
				 x)
			     (x-component (velocity body))))))
    (set-time body time)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Gravitational acceleration on earth = -9.81 m/s^2
(defconstant g (make-vector 0 (m// -981 100)))

;; A 10 kg ball is fired out of a cannon from (0, 0) at 20 m/s at 20
;; degrees above horizontal. When will it hit the ground?
(defparameter *ball*
  (make-instance 'free-body
		 :acceleration g
		 :velocity (make-vector-polar
			    20 (degrees->radians 20))
		 :mass 10))

(coordinates *ball*)
(set-time *ball* 0.5)
(set-time *ball* '$t)
(set-time *ball* #$ 4000*sin(%pi/9)/981 $)

(progn
  (set-x *ball* 2)
  (coordinates *ball*))

;; x(t) := (x-component v)*t + xi
;; y(t) := (1/2)*g*t^2 + (y-component v)*t + yi